function [k1_target_set,k2_target_set,x1_target_set,x2_target_set,grid_x1_target_set,grid_x2_target_set,boo1,boo2,...
    gu_w, gu_target_set,gu_sample_set,dr,dr_hh,Xu_target_set,nw_target_set,n_target_set,Xu_sample_set,nw_sample_set,n_sample_set,x1_sample_set,x2_sample_set,...
    scaled_att_target_set,scaled_att_sample_set,scaled_att_sample_set_weighted,Ind_target_set,Ind_sample_set] ...
    = generate_inputs_cross_validation_quadrant(reweight,sample,sample_training,sample_testing,sample_wave,target_wave,tcost,covariate,SMC)
% This function generates the grid search inputs for the 
% out-of-sample performance exercise 

% Inputs: 
% (1) reweight: use density ratio to reweight waves or not
% (2) sample: cross-sectional dataset for the pooled sample
% (3) sample_training: training sample
% (4) sample_testing: testing sample
% (5) sample_wave: set to 0 to indicate full estimation sample
% (6) target_wave: set to 0 to indicate full estimation sample
% (7) tcost: private marginal cost of implementing the program per household
%  >0 represents cost savings; =0 represents kWh reduction 
% (8) covariate: indicates covariate for analysis: income, size, vintage,
% minimum of baseline consumption, maximum of baseline consumption, or
% standard deviation of consumption
% (9) SMC: =1 use social marginal cost; =0 use retail electricity price

% Outputs: 
% The list is too long to include here. The variables are equivalent to the
% ones used in the main grid search script. 

%% Generate inputs
    wave=0;
    [X,Y,D,n,ps]=generate_inputs_quadrant_rule(sample,wave,tcost,covariate,SMC);
    
if sample_wave==0 && target_wave==0    
    [X_t,Y_t,D_t,n_t,ps_t]=generate_inputs_quadrant_rule_cv(sample_training,wave,tcost,covariate,SMC); 
    [X_h,Y_h,D_h,n_h,ps_h]=generate_inputs_quadrant_rule_cv(sample_testing,wave,tcost,covariate,SMC);
    scaled_att_h= (sum(D_h.*Y_h-((1-D_h).*Y_h.*ps_h)./(1-ps_h))/sum(D_h))*sum(D_h)/n_h;
    scaled_att_t= (sum(D_t.*Y_t-((1-D_t).*Y_t.*ps_t)./(1-ps_t))/sum(D_t))*sum(D_t)/n_t;

elseif sample_wave~=0 && target_wave~=0    

    wave=3;
    [X3_t,Y3_t,D3_t,n3_t,ps3_t]=generate_inputs_quadrant_rule_cv(sample_training{1},wave,tcost,covariate,SMC); 
    [X3_h,Y3_h,D3_h,n3_h,ps3_h]=generate_inputs_quadrant_rule_cv(sample_testing{1},wave,tcost,covariate,SMC);
    scaled_att_wave3_h= (sum(D3_h.*Y3_h-((1-D3_h).*Y3_h.*ps3_h)./(1-ps3_h))/sum(D3_h))*sum(D3_h)/n3_h;
    scaled_att_wave3_t= (sum(D3_t.*Y3_t-((1-D3_t).*Y3_t.*ps3_t)./(1-ps3_t))/sum(D3_t))*sum(D3_t)/n3_t;

    wave=6;
    [X6_t,Y6_t,D6_t,n6_t,ps6_t]=generate_inputs_quadrant_rule_cv(sample_training{2},wave,tcost,covariate,SMC);
    [X6_h,Y6_h,D6_h,n6_h,ps6_h]=generate_inputs_quadrant_rule_cv(sample_testing{2},wave,tcost,covariate,SMC);
    scaled_att_wave6_h= (sum(D6_h.*Y6_h-((1-D6_h).*Y6_h.*ps6_h)./(1-ps6_h))/sum(D6_h))*sum(D6_h)/n6_h;
    scaled_att_wave6_t= (sum(D6_t.*Y6_t-((1-D6_t).*Y6_t.*ps6_t)./(1-ps6_t))/sum(D6_t))*sum(D6_t)/n6_t;

    wave=7;
    [X7_t,Y7_t,D7_t,n7_t,ps7_t]=generate_inputs_quadrant_rule_cv(sample_training{3},wave,tcost,covariate,SMC);
    [X7_h,Y7_h,D7_h,n7_h,ps7_h]=generate_inputs_quadrant_rule_cv(sample_testing{3},wave,tcost,covariate,SMC);
    scaled_att_wave7_h= (sum(D7_h.*Y7_h-((1-D7_h).*Y7_h.*ps7_h)./(1-ps7_h))/sum(D7_h))*sum(D7_h)/n7_h;
    scaled_att_wave7_t= (sum(D7_t.*Y7_t-((1-D7_t).*Y7_t.*ps7_t)./(1-ps7_t))/sum(D7_t))*sum(D7_t)/n7_t;
end
    %% Estimate g using IPW
    g = Y.*(D./ps - (1-D)./(1-ps)); % elements of W(G) - W(G_0)

    if sample_wave==0 && target_wave==0 % pooled 
         g_t = Y_t.*(D_t./ps_t - (1-D_t)./(1-ps_t)); % elements of W(G) - W(G_0), training
         g_h = Y_h.*(D_h./ps_h - (1-D_h)./(1-ps_h)); % elements of W(G) - W(G_0), testing
    elseif sample_wave~=0 && target_wave~=0    
        g3_t = Y3_t.*(D3_t./ps3_t - (1-D3_t)./(1-ps3_t)); 
        g6_t = Y6_t.*(D6_t./ps6_t - (1-D6_t)./(1-ps6_t)); 
        g7_t = Y7_t.*(D7_t./ps7_t - (1-D7_t)./(1-ps7_t)); 
        g3_h = Y3_h.*(D3_h./ps3_h - (1-D3_h)./(1-ps3_h)); 
        g6_h = Y6_h.*(D6_h./ps6_h - (1-D6_h)./(1-ps6_h)); 
        g7_h = Y7_h.*(D7_h./ps7_h - (1-D7_h)./(1-ps7_h)); 
    end
    %% Re-weight waves by wave density ratio

    % sort rows of X -> [Xr, gr]
    [Xr, Ind] = sortrows(X);
    gr = g(Ind);
    
    if sample_wave~=0 && target_wave~=0
        [Xr3_h, Ind3_h] = sortrows(X3_h);
        gr3_h = g3_h(Ind3_h);
        [Xr6_h, Ind6_h] = sortrows(X6_h);
        gr6_h = g6_h(Ind6_h);
        [Xr7_h, Ind7_h] = sortrows(X7_h);
        gr7_h = g7_h(Ind7_h);
    elseif sample_wave==0 && target_wave==0
        [Xr_h, Ind_h] = sortrows(X_h);
        gr_h = g_h(Ind_h);
    end
        
    switch covariate
        case 'income'
            N_bins_x=29; 
            N_bins_y=272; 
        case 'size'
            N_bins_x=48; 
            N_bins_y=272; 
        case 'vintage'
            N_bins_x=17; 
            N_bins_y=272; 
    end

    % generate grid lower boundusing the full sample
    [N,boo1,boo2,xx1,xx2]=histcounts2(Xr(:,1),Xr(:,2),...
        [N_bins_x,N_bins_y],'Normalization','count'); 
    %boo1: the lower bound of each bin for income
    %boo2: the lower bound of each bin for prev_use
    %xx1: each data point belongs to which income bin 
    %xx2: each data point belongs to which prev_use bin 

    nu = length(boo1)*length(boo2);
    % number of unique grid point= number of income bin * number of prev_use

    if sample_wave~=0 && target_wave~=0
        [N_3,boo1_3,boo2_3,xx1_3,xx2_3]=histcounts2(Xr3_h(:,1),Xr3_h(:,2),...
            'XBinEdges',boo1,'YBinEdges',boo2);
        [N_6,boo1_6,boo2_6,xx1_6,xx2_6]=histcounts2(Xr6_h(:,1),Xr6_h(:,2),...
            'XBinEdges',boo1,'YBinEdges',boo2);
        [N_7,boo1_7,boo2_7,xx1_7,xx2_7]=histcounts2(Xr7_h(:,1),Xr7_h(:,2),...
            'XBinEdges',boo1,'YBinEdges',boo2); 

        [Xu3_h,gu3_h,nw3,x1_wave3,x2_wave3] = empirical_density_generate_grid(nu,boo1,boo2,xx1_3,xx2_3,gr3_h);
        [Xu6_h,gu6_h,nw6,x1_wave6,x2_wave6] = empirical_density_generate_grid(nu,boo1,boo2,xx1_6,xx2_6,gr6_h);
        [Xu7_h,gu7_h,nw7,x1_wave7,x2_wave7] = empirical_density_generate_grid(nu,boo1,boo2,xx1_7,xx2_7,gr7_h);

        grid_x1_wave3=boo1;
        grid_x1_wave6=boo1;
        grid_x1_wave7=boo1;
        grid_x2_wave3=boo2;
        grid_x2_wave6=boo2;
        grid_x2_wave7=boo2;

        k1_wave3=N_bins_x+1;
        k1_wave6=N_bins_x+1;
        k1_wave7=N_bins_x+1;

        k2_wave3=N_bins_y+1;
        k2_wave6=N_bins_y+1;
        k2_wave7=N_bins_y+1;
        
    elseif sample_wave==0 && target_wave==0
        [N_h,boo1_h,boo2_h,xx1_h,xx2_h]=histcounts2(Xr_h(:,1),Xr_h(:,2),...
            'XBinEdges',boo1,'YBinEdges',boo2);
        [Xu_h,gu_h,nw_h,x1_h,x2_h] = empirical_density_generate_grid(nu,boo1,boo2,xx1_h,xx2_h,gr_h);
        grid_x1_h=boo1;
        grid_x2_h=boo2;
        k1_h=N_bins_x+1;
        k2_h=N_bins_y+1;
    end
    % Define x1,x2,grid_x1,grid_x2 for target population, which will be used in
    % the grid search part 
    if target_wave==3
        x1_target_set=x1_wave3;
        x2_target_set=x2_wave3;
        grid_x1_target_set=grid_x1_wave3;
        grid_x2_target_set=grid_x2_wave3;
        k1_target_set=k1_wave3;
        k2_target_set=k2_wave3;
        Xu_target_set=Xu3_h;
        nw_target_set=nw3;
        n_target_set=n3_h;
        gu_target_set=gu3_h;
        scaled_att_target_set=scaled_att_wave3_h;
        Ind_target_set=Ind3_h;

    elseif target_wave==6
        x1_target_set=x1_wave6;
        x2_target_set=x2_wave6;
        grid_x1_target_set=grid_x1_wave6;
        grid_x2_target_set=grid_x2_wave6;
        k1_target_set=k1_wave6;
        k2_target_set=k2_wave6;
        Xu_target_set=Xu6_h;
        nw_target_set=nw6;
        n_target_set=n6_h;
        gu_target_set=gu6_h;
        scaled_att_target_set=scaled_att_wave6_h;
        Ind_target_set=Ind6_h;


    elseif target_wave==7
        x1_target_set=x1_wave7;
        x2_target_set=x2_wave7;
        grid_x1_target_set=grid_x1_wave7;
        grid_x2_target_set=grid_x2_wave7;
        k1_target_set=k1_wave7;
        k2_target_set=k2_wave7;
        Xu_target_set=Xu7_h;
        nw_target_set=nw7;
        n_target_set=n7_h;
        gu_target_set=gu7_h;
        scaled_att_target_set=scaled_att_wave7_h;
        Ind_target_set=Ind7_h;
        
    elseif target_wave==0
        x1_target_set=x1_h;
        x2_target_set=x2_h;
        grid_x1_target_set=grid_x1_h;
        grid_x2_target_set=grid_x2_h;
        k1_target_set=k1_h;
        k2_target_set=k2_h;
        Xu_target_set=Xu_h;
        nw_target_set=nw_h;
        n_target_set=n_h;
        gu_target_set=gu_h;
        scaled_att_target_set=scaled_att_h;
        Ind_target_set=Ind_h;
    end

    %% Define gu for sample wave 
    if sample_wave~=0 && target_wave~=0
        [Xr3_t, Ind3_t] = sortrows(X3_t);
        gr3_t = g3_t(Ind3_t);
        [Xr6_t, Ind6_t] = sortrows(X6_t);
        gr6_t = g6_t(Ind6_t);
        [Xr7_t, Ind7_t] = sortrows(X7_t);
        gr7_t = g7_t(Ind7_t);

        [N_3,boo1_3,boo2_3,xx1_3,xx2_3]=histcounts2(Xr3_t(:,1),Xr3_t(:,2),...
            'XBinEdges',boo1,'YBinEdges',boo2);
        [N_6,boo1_6,boo2_6,xx1_6,xx2_6]=histcounts2(Xr6_t(:,1),Xr6_t(:,2),...
            'XBinEdges',boo1,'YBinEdges',boo2);
        [N_7,boo1_7,boo2_7,xx1_7,xx2_7]=histcounts2(Xr7_t(:,1),Xr7_t(:,2),...
            'XBinEdges',boo1,'YBinEdges',boo2); 

        [Xu3,gu3,nw3,x1_wave3,x2_wave3] = empirical_density_generate_grid(nu,boo1,boo2,xx1_3,xx2_3,gr3_t);
        [Xu6,gu6,nw6,x1_wave6,x2_wave6] = empirical_density_generate_grid(nu,boo1,boo2,xx1_6,xx2_6,gr6_t);
        [Xu7,gu7,nw7,x1_wave7,x2_wave7] = empirical_density_generate_grid(nu,boo1,boo2,xx1_7,xx2_7,gr7_t);
    
    elseif sample_wave==0 && target_wave==0
        [Xr_t, Ind_t] = sortrows(X_t);
        gr_t = g_t(Ind_t);
        [N_3,boo1_t,boo2_t,xx1_t,xx2_t]=histcounts2(Xr_t(:,1),Xr_t(:,2),...
            'XBinEdges',boo1,'YBinEdges',boo2);
        [Xu_t,gu_t,nw_t,x1_t,x2_t] = empirical_density_generate_grid(nu,boo1,boo2,xx1_t,xx2_t,gr_t);
    end
    
    % get sample wave parameters for export
    if sample_wave==3
        x1_sample_set=x1_wave3;
        x2_sample_set=x2_wave3;
        Xu_sample_set=Xu3;
        nw_sample_set=nw3;
        n_sample_set=n3_t;
        gu_sample_set=gu3;
        scaled_att_sample_set=scaled_att_wave3_t;
        Ind_sample_set=Ind3_t;

    elseif sample_wave==6
        x1_sample_set=x1_wave6;
        x2_sample_set=x2_wave6;
        Xu_sample_set=Xu6;
        nw_sample_set=nw6;
        n_sample_set=n6_t;
        gu_sample_set=gu6;
        scaled_att_sample_set=scaled_att_wave6_t;
        Ind_sample_set=Ind6_t;

    elseif sample_wave==7
        x1_sample_set=x1_wave7;
        x2_sample_set=x2_wave7;
        Xu_sample_set=Xu7;
        nw_sample_set=nw7;
        n_sample_set=n7_t;
        gu_sample_set=gu7;
        scaled_att_sample_set=scaled_att_wave7_t;
        Ind_sample_set=Ind7_t;
        
    elseif sample_wave==0
        x1_sample_set=x1_t;
        x2_sample_set=x2_t;
        Xu_sample_set=Xu_t;
        nw_sample_set=nw_t;
        n_sample_set=n_t;
        gu_sample_set=gu_t;
        scaled_att_sample_set=scaled_att_t;
        Ind_sample_set=Ind_t;
    end
    
    %% calculate density ratio and weighted gu 
    if reweight==0
        gu_w=gu_sample_set;
        dr_hh=ones(n_t,1); 
        dr=ones(size(nw_t,1),1);
        scaled_att_sample_set_weighted= scaled_att_sample_set;

    else
        dr=((nw_target_set)./n_target_set)./(nw_sample_set./n_sample_set); 
        dr(isinf(dr))=0;
        dr(isnan(dr))=0;
        gu_w=gu_sample_set.*dr;

    % get household specific density ratio, then calculate
    % scaled_att_sample_set weighted by density ratio
    if sample_wave==3
        dr_hh = nan(size(g3_t,1),1); 
        i_u=1;

            for i=1:length(boo1_3)
                for j=1:length(boo2_3)
                    bin_boo = logical((xx1_3==i).*(xx2_3==j));
                    dr_ = dr(i_u);
                    dr_hh(bin_boo,1) = dr_;
                    i_u=i_u+1;
                end
            end
        scaled_att_sample_set_weighted= (sum(D3_t(Ind3_t).*Y3_t(Ind3_t).*dr_hh-((1-D3_t(Ind3_t)).*Y3_t(Ind3_t).*dr_hh.*ps3_t)./(1-ps3_t))/sum(D3_t(Ind3_t)))*sum(D3_t(Ind3_t))/n3_t; % Ind is the sorted order in the pooled sample
    elseif sample_wave==6
        dr_hh = nan(size(g6_t,1),1); 
        i_u=1;

            for i=1:length(boo1_6)
                for j=1:length(boo2_6)
                    bin_boo = logical((xx1_6==i).*(xx2_6==j));
                    dr_ = dr(i_u);
                    dr_hh(bin_boo,1) = dr_;
                    i_u=i_u+1;
                end
            end
        scaled_att_sample_set_weighted= (sum(D6_t(Ind6_t).*Y6_t(Ind6_t).*dr_hh-((1-D6_t(Ind6_t)).*Y6_t(Ind6_t).*dr_hh.*ps6_t)./(1-ps6_t))/sum(D6_t(Ind6_t)))*sum(D6_t(Ind6_t))/n6_t; % Ind is the sorted order in the pooled sample
    elseif sample_wave==7
        dr_hh = nan(size(g7_t,1),1); 
        i_u=1;

            for i=1:length(boo1_7)
                for j=1:length(boo2_7)
                    bin_boo = logical((xx1_7==i).*(xx2_7==j));
                    dr_ = dr(i_u);
                    dr_hh(bin_boo,1) = dr_;
                    i_u=i_u+1;
                end
            end
        scaled_att_sample_set_weighted= (sum(D7_t(Ind7_t).*Y7_t(Ind7_t).*dr_hh-((1-D7_t(Ind7_t)).*Y7_t(Ind7_t).*dr_hh.*ps7_t)./(1-ps7_t))/sum(D7_t(Ind7_t)))*sum(D7_t(Ind7_t))/n7_t; % Ind is the sorted order in the pooled sample
    end
    end
end